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We propose fast, exact and efficient algorithms for the convolution of two arbitrary functions on 
the sphere which speed up computations by a factor O^y/Wj compared to present methods where 
N is the number of pixels. No simplifying assumptions are made other than bandlimitation. This 
reduces typical computation times for convolving the full sky with the asymmetric beam pattern of 
a megapixel Cosmic Microwave Background (CMB) mission from months to minutes. Our methods 
enable realistic simulation and careful analysis of data from such missions, taking into account the 
effects of asymmetric "point spread functions" and far side lobes of the physical beam. While 
motivated by CMB studies, our methods are general and hence applicable to the convolution or 
filtering of any scalar field on the sphere with an arbitrary, asymmetric kernel. We show in an 
appendix that the same ideas can be applied to the inverse problems of map-making and beam 
reconstruction by similarly accelerating the transpose convolution which is needed for the iterative 
solution of the normal equations. 



I. INTRODUCTION 

A major near-term objective in the field of Cosmol- 
ogy today is to gain a detailed measurement and statis- 
tical understanding of the anisotropics of the cosmic mi- 
crowave background (CMB). While the theory of primary 
CMB anisotropy is well-developed (see Jl|] for a review) 
and we are facing a veritable flood of data from a new 
generation of instruments and missions, perhaps the sin- 
gle most limiting factor for interpreting these data is the 
exorbitant computational cost involved in realistic mis- 
sion simulation and careful analysis of the data products 

II- 

Important and computationally expensive tasks for 
both simulation and analysis of microwave data are to 
simulate and to correct for the systematic errors due 
to imperfections of realistic microwave telescopes, such 
as beam asymmetries and far side lobes. The effect of 
an asymmetric "point spread function" is to distort the 
shapes of the detected anisotropies. What makes far 
side lobes an important issue is the fact that the CMB 
anisotropy signal has an amplitude of one in 10 5 rela- 
tive to the 2.7 K background. In regions of low galac- 
tic latitude, foregrounds from galactic synchrotron radi- 
ation and dust emission are expected to exceed this signal 
by many orders of magnitude over a wide range of fre- 
quencies . Even though CMB experiments will obvi- 
ously not target these regions to obtain measurements of 
the background anisotropy, the large amplitudes of these 
galactic sources may induce systematic errors even when 
"looking" in directions far away from the galactic plane 
if the instrument allows diffraction of stray light into the 
detectors. Solar system bodies, including the earth, are 
other possible sources of stray light. 



To assess these problems and formulate solutions we 
must be able to compute the detector response at ev- 
ery pointing of the telescope. The inputs are a physi- 
cal model of the "beam" over An steradian and a model 
of the "sky" containing both simulated signal as well as 
foreground sources possibly including ground emission. 
Note that in the general case not just the direction of 
the pointing is important but also the orientation of the 
beam about the pointing axis. The detector response is 
then the solution to a quadrature problem at each orien- 
tation. 

Analysis methods of CMB data have neglected this 
difficulty by assuming azimuthal symmetry of the beam 
which greatly simplifies the calculation [^|-^| . Simulation 
work which did include an asymmetric beam and far side 
lobes using pixel based methods J|-[ll| ran up against 
computational challenges for angular scales smaller than 
one degree, running for hundreds of hours even with 
optimized adaptive mesh algorithms. Such algorithms 
are clearly inadequate for modern high resolution exper- 
iments which achieve resolutions of a few minutes of arc. 

In this paper we describe a numerical method which 
greatly accelerates the computations which are necessary 
to correctly account for realistic beam profiles in simula- 
tion and analysis of directional data on the sphere. This 
is achieved by rewriting the problem in such a way that 
we can take advantage of the Cooley-Tukey Fast Fourier 
Transform (FFT) algorithm. 

The following section of this paper defines the general 
problem in terms of rotations of the beam with respect 
to the sky. We then introduce a geometrically motivated 
split of the rotation operator in section three. This en- 
ables us, in section four, to derive the general solution 
for the detector response for all possible relative orienta- 
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tions of the beam and the sky within a given section on 
the sphere. Section five then discusses the solution and 
derives special cases from it, amongst others the well- 
known algorithm for convolution with azimuthally sym- 
metric kernels. We conclude in section six with a timing 
example. An appendix applies the same ideas to acceler- 
ating the computation of the transverse convolution, an 
operation which becomes important in the inverse prob- 
lem of map estimation. 

While we are motivated by the goal of achieving and 
interpreting precision measurements of the anisotropics 
of the cosmic microwave background, the methods wc 
present are general and apply to the convolution or fil- 
tering of any scalar field on the sphere with an arbitrary, 
asymmetric but constant kernel. We generalise our meth- 
ods to tensor fields on the sphere in reference fll2[]. 



simple operation in position space (the pixel basis) and 
can be computed separately. Linearity allows us to then 
add the results of the m convolution of extended sources 
to the point source convolution. 

In the most general case, the bandlimits (see Eq. (^) 
below for a definition) of the beam and the sky are 
and L s , respectively. Define L = min(Lf,,L s ). Note that 
we actually only need one of s, b to be bandlimited as 
long as the multipoles of the other are bounded as / — > 
oo. Then the numerical evaluation of the integral in Eq. 
([[]) takes C(L 2 ) operations for each tuple ($1, 6, <& 2 ). 
These integrals need to be evaluated for a grid of beam 
locations that has to contain at least C?(L 3 ) grid points 
to allow subsequent interpolation at arbitrary locations. 
Therefore the total computational cost for evaluating the 
convolution using Eq. (Q) scales as C(L 5 ) . 



II. STATEMENT OF THE PROBLEM 



III. FACTORIZING THE ROTATION 



Consider two bandlimited functions on the sphere 6(7) 
and 5(7). For definiteness and to aid the imagination we 
will refer to them in the following as the beam and the 
sky, respectively, but they could be completely general 
bandlimited functions in particular neither of them is 
constrained to be positive definite or even real. 

The task is to compute the scalar product of the beam 
and the sky at a set of beam orientations. To describe 
these orientations, we use the Euler angles <E>i,0 and 
$20. The convolved signal for each beam orientation 
($1, O, §2) can then be written as 

T($ 2 ,e,$!) = J <% [5(# 2j e, $i)b\(i)*s(i). (i) 

Here the integration is over all solid angles, D is the 
operator of finite rotations such that Db is the rotated 
beam, and the asterisk denotes complex conjugation. 

If ($1, 0, $2) can be written as a continuous function 
of a parameter t S [0, T], say, then we call the ordered 
set of tuples ($i(t), 0(i), ( i> 2 (<)) a scan path. Note that 
Eq. |TJ) assumes that time varying signals in the sky vary 
either on time scales much longer than the duration of 
the scan or much smaller than the integration time per 
sample. In the context of CMB missions this is a good 
approximation with the exceptions of planets (for long 
duration missions), time varying point sources, and at- 
mospheric foregrounds. Of these only atmospheric fore- 
grounds present a problem for the convolution, because 
they are extended - convolution with a point source is a 



It is possible to simplify the evaluation of Eq. ([j]) sig- 
nificantly by factorizing the rotation into two auxiliary 
rotations such as 



*Our Euler angle convention refers to active right handed 
rotations of a physical body in a fixed coordinate system. 
The coordinate axes stay in place under all rotations and the 
object rotates around the z, y and z sixes by $1, O and $2, 
respectively, according to the right handed screw rule. 



D($ 2 , 0, $1) = D(</> E , E , 0)D((f>, 8, to). 



(2) 



We will define the various angles and motivate this split 
in the following. Figure [l] is intended to illustrate this 
discussion. 

To introduce these coordinates let us first consider ba- 
sic scan paths. Imagine a scan path where the beam 
sweeps over the sky by scanning on rings of angular ra- 
dius 8 £ [0, 7r/2). The centers of these scanning circles lie 
on a ring of constant latitude, a polar angle 6e € [0, tt) 
away from the north pole. The angle cf>E G [0, 2-tt) selects 
a given scanning circle and is defined as the longitude of 
its center, while <f> G [0, 27r) measures the angle along each 
scanning ring defined as increasing in a right-handed way 
about the outward normal at the center, starting from 
zero at the southernmost point on the ring. Hence, for 
such a basic scan path we can write the convolution as 
set of scalar products T (</>£, </>). The angles 8 and 6e 
are thought of as parameters which fixed to define the 
scanning geometry. 

As a generalization of basic scan paths, we allow as 
a further degree of freedom an additional right handed 
rotation of the beam about its outward axis by an angle 

lu e [0,7r/2). 

Now we can see geometrically that all beam orienta- 
tions on generalized basic scan paths can be arrived at by 
successively applying the two rotations in Eq. (|J). Define 
as the null position of the beam when it is oriented along 
the z-axis 8 — 8e — 4>e — 4> — u — 0. Starting from 
the null position, acting on it with D((j),8,uj) rotates it 
about its axis by ui and moves it out onto a ring with 
opening angle 8 at an azimuthal angle <fi. Then acting 
with D((I)e,8e,0) moves the beam into position. 
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l=L f m =l 

/Of) = Y flmYlmil), 

1=0 m=-l 



(4) 



where 7 denotes a unit vector. For practical applica- 
tions the bandlimit Lf is chosen such that higher terms 
contribute insignificantly. We use the notation where all 
quantities carrying both an I and an m index vanish for 
m > I. This saves having to write explicit limits for sums 
over azimuthal quantum numbers. 

Invariance of the scalar product under a change of basis 
then allows us to re-write Eq. (pi) as 



T(<f> E) <i>,u)) = 



I III 



D( 



Q)D{(j>,e,w)b 



! 111 



(5) 



Y s hnD l n * M ((f>E, E , 0)D 1 m M , 0, 9, uj)b* M ,. 



ImMM' 



FIG. 1. Our coordinate system for efficient convolutions. 
The beam is shown at the position corresponding to 9 = 35°, 
6e — 50°, (f>E = 60°, 4> — 0° and u> = 0°. The cross-hairs on 
the beam mark its orientation, here shown for uo = 0. In the 
null position (8 = 9e = <t>E = <t> = w — 0) the beam is aligned 
with the z-axis, the vertical cross hair pointing along increas- 
ing x and the horizontal cross hair pointing along increasing 

y- 



Using the factorization Eq. (Eh we can re- write Eq. (hi) 



T(cf > E,<j>,uj)= dn=f D(</>E,0E,O)D(<l>,e,u)b(i)*s(i), 



A simple explicit expression for the matrix elements 
D l mm ,(cf>2,0, 4>i) can be given. One can define a real func- 
tion d l mm , (9) such that 



(6) 



Thus the dependence of D on the Euler angles <pi and <p2 
is only in terms of complex exponentials. While explicit 
formulas for the d- functions exist |l3|| , they are more 
conveniently computed numerically using their recursion 
properties [[LjJ . 

Substituting into Eq. (|J) and defining the three- 
dimensional Fourier transform of T(<J)e, as 



9, B E fixed). (3) ( 2 7r) 2 



d(j) E d(j)dw T(4> E , <t>, u) e -*m*»-*»V- 



The function T{4>e, 4>, oj) contains all possible integrals 
for a given scanning geometry. In fact, for the special 
case 9 = 9e = vr/2, these angles parameterize all pos- 
sible orientations of the beam on the sky, i.e. (<^e,(/>, oj) 
parameterize the group of rotations in three dimensions. 
It is well known that in this case these coordinates cover 
SO (3) twice, but this can be easily remedied by restrict- 
ing the range of one of the angles to half its range. We 
defer removing this redundancy until the end of our cal- 
culation. 



IV. SOLUTION 

To exploit the form of Eq. (g), it is expedient to repre- 
sent the functions s and 6 as well as the rotation operators 
in the spherical harmonic basis. A bandlimited function 
f(ff) can be expanded in spherical harmonics as 



we obtain 



T 



m m m 



I 



Sir, 



iMi^E)d MM ,(0)b* M , 



(7) 



(8) 



This equation is the main result of this paper, in ef- 
fect generalising fast 2D Fourier transfrom convolution 
from the plane to the sphere. Its properties and spe- 
cialisations will be discussed in the next section. Here 
we give a geometrical interpretation. We have arrived at 
Eq. (^) by writing convolution problems in such a way 
that the results are fields on 3-tori instead of subsets of 
the 3-sphere, which is the group manifold of rotations in 
3 dimensions. Convolutions over azimuthally symmetric 
and connected sections of the 2-sphere (such as polar caps 
or annuli) can be parameterised by 9 and 9e and hence 
can be extended to 3-tori as shown. Since exponentials 
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are a complete and orthonormal basis on the 3-torus and 
because we assumed that s and b are band-limited, the 
T m m ' m" contain all information about the inverse trans- 
form, Eq. |), 



T((f> E ,<j>,Uj) 



m,m',m"--L 



imcpE -Vim! <j}-Vim" to 



(9) 



Not all tuples (4>e, <j>, uj) correspond to distinct beam ori- 
entations but this redundancy is more than compensated 
for by the efficiency of the method. 



V. DISCUSSION 



Several remarks about Eq. (^|) are in order. 



1. Computational cost 



Computing the T m m i m ii in Eq. 

(|J) costs C(L 4 | svo.9e\29 /tt) operations. The factors 
9 and |sin^| come from the fact that the bandlimit L 
for the sky implies a bandlimit oc 29L/ir on a ring of ra- 
dius 9 and hence the ranges of m and m' can be reduced 
by factors of | sin 9e\ and 29 /it, respectively, if the rings 
and the ring of ring centers are not great circles. Using 
the Fast Fourier Transform (FFT) algorithm, the inverse 
Fourier transform takes O ( L 3 log L | sin 9e \ 20 /ti) oper- 
ations. If the convolved sky is assumed to be real, we 
have 



- m m' m" 



— m — 7n' —m" ! 



(10) 



reducing memory and processor requirements by a factor 
2. 



2. Quadrature and interpolation 

In pixel space each evaluation of T(cf>E,4>, uj) is an ex- 
plicit quadrature problem and hence necessarily approx- 
imate. In our approach, all sums have a finite number 
of terms and the results are exact as long as s and b are 
band-limited. Quadrature issues only have to be dealt 
with if 6 or s are given in pixel space and we have to 
evaluate the beam and sky multipole coefficients bi m and 
si m . The details of which pixelization to choose on the 
sphere and how to solve this generalized quadrature prob- 
lem for the multipole coefficients are outside of the scope 
of this work but an efficient and practical approach to 
the quadrature problem is implemented in the HEALPix 
package and will be discussed in a future publication 



An interesting property of Eq. (||) is that as long as L 
was chosen appropriately one is guaranteed to have the 
convolved sky sampled sufficiently densely for worry-free 
interpolation on either of the three indices. 



A. Special cases 
We will now discuss certain special cases of Eq. (^) . 

1. Total convolution 

Let us obtain the convolved sky at all possible beam 
orientations to on an equidistant coordinate grid in 4> (cor- 
responding to the polar angle) and <pE (corresponding to 
the azimuthal angle) . We will refer to this case as the to- 
tal convolution. This can be achieved by evaluating Eq. 
(||) setting 9 = 6e = f ■ In this case we only need to 
know dj„/ m (f )■ This means we only have to evaluate a 
single recursion relation to evaluate the sum on I, which 
simplifies the computation. The inverse FFT gives the 
desired result. 

A further simplification arises in this case from the fact 
that if 9e = § , not all components of T m m / m i> are inde- 
pendent. The redundancy in the parametrization where 
the polar angle 4> € [0, 2ir) leads to the symmetry 



Tt 



j) — T(tt + <j) E , 2n - (f), tt + uj). 



(11) 



This translates into the identity 



± mm ! m" — V x J J -m—'m'm"i \ ±Zj J 

which cuts the required memory and computation time 
by a factor 2. 



2. Exact or approximate azimuthal symmetry of the beam 

In many practical situations the "beam" represents the 
response function of an optical system with only mild im- 
perfections. If this is the case, the beam has only slowly 
varying azimuthal structure, implying a cutoff wavenum- 
ber M such that bi m ~ for m > M. in this case 
the computational cost for a total convolution scales as 
O^Mflsin^) . 

In the limit of an azimuthally symmetric beam, M = 0, 
we obtain an 0(L 3 9sm9E) method. However, it is 
known [jl5| that at least in principle there exist faster 
methods for convolution of a function on the full sphere 
{9 = 9e — tt/2) with an azimuthally symmetric beam 
which scale as 0(L 2 (log(LV) 2 ) . We can show how this 
limit is obtained from Eq. (0) by using the facts that the 
^mm'(§ ) are the Fourier coefficients of the d l mm ,(9) and 
that d! mO (0) = Pi m {9). Then Eq. (§) reduces to the form 



0- 



4 



T(<f> E ,tf>) 



, YlmilT ~(j>,4>E+ ir/2)bi sir. 



(13) 



where the arguments of Yi m are the polar angle and 
the azimuthal angle, respectively. The algorithm by |l5| ] 
succeeds precisely in reducing the computational cost of 
evaluating this expression to C(L 2 (log 2 L)) under the 
proviso of the technical difficulties there outlined. 

We note here for completeness, that by choosing a delta 
function beam (and hence 6; TO = const), we recover the 
Fourier summation method for the spherical harmonic 
transform, described in equations (5.2) to (5.4) in [ p^[ . 
This computes Eq. (Q) on an equidistant coordinate grid 
by doing Fourier transforms on latitudinal and longitu- 
dinal lines. The forward transform is obtained by simply 
working all steps in reverse. 



3. Basic scan paths 

Consider an application where the convolution is re- 
quired only along a "basic scan path" . This is one of the 
proposed scan strategies for the Planck satellite mission. 
From our definition of basic scan paths in section III we 
see that they correspond to setting uj 



in Eq. (g). 

Computing the inverse Fourier transform of Eq. (||) 
with u> — just amounts to summing over m". Then 
only the two-dimensional Fourier transform 



,(u> = 0) = 



E 



s lmd mm , 



hn' 



remains to be evaluated. The quantity 
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(14) 



(15) 



M 



can be precomputed. All in all the computational time 
needed for evaluating these expressions is O(L 3 sin Oe)- 
Storage requirements scale only as C(L 2 ) . 

Note that in this case azimuthal symmetry does not 
necessarily imply reduced computational cost. If the 
beam is concentrated at the north pole into a region 
of size a then Xi m will have m modes populated up to 
M~e/a. 

Geometrically, the basic scan path corresponds to a 
2-torus which is the section of the 3-sphere of rotations 
at constant u>. Note that in this case there is no redun- 
dancy in the parametrization — every tuple (4>e, 0,0) 
corresponds to a distinct beam orientation. 



4- Perturbations about basic scan paths 



Such scan paths can be composed by computing sev- 
eral convolutions along basic scan paths for different an- 
gles 9e and then choosing scanning circles at will from 
among these basic ones. This method suggests itself if 
the precession angle is small and hence a small number 
of convolutions is sufficient to sample the variation in Oe- 
Convolutions at points which do not coincide with sam- 
pling points can then be determined by interpolation. 

Another approach to this type of problem and further 
generalisations are discussed in the next paragraphs. 



5. Other special cases 

Other potentially interesting special cases of Eq. (||) 
can be worked out by fixing any of the parameters to spe- 
cial values and evaluating the inverse transform, analo- 
gous to the calculation for basic scan paths. For example 
one obtains expressions for 

• All possible beam orientations along a circle of con- 
stant latitude Oe- In this case <j> and u> have the 
same meaning; formally, T mm / TO » = T mm /<5 m / m /> 
and we obtain an C(L 2 M sin^s) method: 



i 



s lmd mm 



(16) 



• Individual scanning rings of a basic scan path. 
Here, u) = 0, and the only free parameter is <f>. 

The details of the calculations for this and similar cases 
are now easy exercises. 

6. Generalizations 

Further, it is clear from the derivation that more gen- 
eral types of paths can be constructed by factorizing the 
rotation operator more than twice, so as to generate for 
example a ring of ring of rings, etc. For particular appli- 
cations some of these may be advantageous, for example 
if they simplify the interpolation problem on the output 
ring set. A specific example is the precessing scan path 
mentioned in the previous subsection. Inserting another 
rotation operator D(0,0p,(f>p) between the two opera- 
tors in Eq. (^), and setting uj = produces a set of rings 
whose centers lie on circles of radius 9p about thetciE- 
This may simplify the interpolation problem. The rota- 
tion corresponding to <j)p can be sampled sparsely if Op 
is small (Mp ~ sin Op sin^L with obvious notation) and 
the interpolation problem becomes simpler. 



A slight generalization of the previous case are scan 
paths which are close to basic but include a variation 
in Oe from scanning circle to scanning circle. Such paths 
result for example from precessing or "wobbling" the spin 
axis of a scanning satellite. 



VI. CONCLUSIONS 

This paper presents a general algorithm which greatly 
reduces the computational cost of convolving two ban- 
dlimited but otherwise arbitrary functions on the sphere. 
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The speedup increases linearly with the smallest angular 
scale of the smoother of the two functions in the prob- 
lem. The scalings of the necessary operation counts are 
discussed in detail in section five. 

We quote in the appendix formulas showing how the 
ideas presented in this paper can be applied to the inverse 
problem of "deconvolution" by speeding up the iterative 
solution of the normal equation in an analogous way. 

This paper focuses on the convolution of scalar valued 
functions on the sphere such as temperature, elevation, 
etc. In order to be able to deal with the polarization of 
the cosmic microwave background we extend the methods 
presented here to tensor valued functions on the sphere 
in reference |l2|] . 

The algorithms which are presented here are already 
being used as a core component of the prototype sim- 
ulation pipeline of the Planck satellite. To give an ex- 
ample for the timing gains one makes by applying this 
method, we computed the following case: both sky and 
beam were interpolated and pixelised very densely, with 
millions of pixels each to resolve the steep variations over 
many orders of magnitude. The bandlimit was somewhat 
generously chosen as L = 1024. Then the convolution of 
the sky with the beam of a single detector for a whole 
year of mission data, consisting of (2049) 2 ~4x 10 6 con- 
volved samples along a basic scan path was generated in 
less than 15 minutes on a single Silicon Graphics R10000 
processor. This compares with several days of computa- 
tion on a severely coarsened sampling grid with several 
hundred times fewer samples on the same machine, using 
the adaptive mesh method Q. For the same resolution 
which we achieved with our methods, the adaptive mesh 
code would have run for months. 

Due to our methods, future CMB missions can go be- 
yond having to approximate the treatment of realistic 
beams. Our methods lend themselves to being used in 
conjunction with iterative map-making methods to re- 
move from the data artefacts which are due to beam dis- 
tortions and far side lobes (see appendix). 

Lastly, we feel that the geometric constructions, analo- 
gies to group properties and algebraic results we intro- 
duce in this article may be useful more generally for CMB 
data analysis and plan to explore these issues in future 
work. 



APPENDIX A: THE CONVOLUTION 
TRANSPOSE 

We sketch here how to set up the inverse problem of 
reconstructing the true sky from the convolved observa- 
tions. If we start with a noise free set of convolutions, 
then the equation to be inverted in order to estimate the 
true underlying sky is, schematically, 



As 



(Al) 



Here, s is the true sky, d is the vector containing the time- 
ordered data after observation. The convolution operator 
is represented by A. 

The least-square estimator for the true sky, § then sat- 
isfies the normal equation 



A As = A d . 



(A2) 



For a perfect observation with a delta function beam, 
A T A = 1. So it may be reasonable to expect that we 
can make progress by considering a mildly imperfect op- 
tical system and consider iterative techniques for solving 
the normal equation. In this case the ability to solve for 
s iteratively (e.g. using a Conjugate Gradient technique) 
relies on convergence (which is assured up to numerical 
effects because the normal matrix A T A is positive def- 
inite and being able to compute the matrix products in 
( |A2j ) quickly. The application of A can be computed ef- 
ficiently using the formulae set out in sections four and 
five. We now find an algorithm for the efficient applica- 
tion of A T , the transpose convolution. 



1. Applying the transpose convolution 

We can write down the expression for A T in a similar 
way to Eq. (||) 



t>Ed<pduj 



D(<l>E,6E,0)D(<t>, e, io)bp)* T(<t> E , 0, w). 

(A3) 
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Now the derivation is analogous to the one preceding Eq. 
(H), and 2/(7) is given in terms of the Wigner d functions 
as 



2/(7)= E Y iW)dL„AeE)d l mlmll (0EK n ,,T, 



mm' in" 



Imm'm" 



(A4) 

This formula can be generalised or applied to special 
cases just as we showed in section |y| for Eq. (||) . 
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